function result = GetBeamPropagation(obj, x, y, z, rateRelax, varargin)
    p = inputParser;
    p.addParameter('output', 'all', @(x) ismember(x, {'all', 'phi', 'rate', 'intensity'}));
    p.parse(varargin{:});
    
    pumpBeam = obj.getPumpBeam();
    absorption = obj.GetCellAbsorption('Pump', 'frequencyRange', [], 'isPlot', false);
    
    cs = absorption.cross_section_average; % [cm^2]
    od = absorption.optical_depth_average; 
    phi0 = pumpBeam.photonFlux(x, y); % [1/s/cm^2]
    alpha0 = phi0*cs/rateRelax;
    
    res.phiZ = (rateRelax/cs).*lambertw(alpha0.*exp(alpha0-od.*z(:)./obj.conditions.cellSize));
    res.rateZ = cs.*res.phiZ;
    res.intensityZ = res.phiZ*(2*pi*h_bar*pumpBeam.parameter.frequency);

    switch p.Results.output
        case 'all'
            result = res;
        case 'phi'
            result = res.phiZ;
        case 'rate'
            result = res.rateZ;
        case 'intensity'
            result = res.intensityZ;
    end
        
end

